Bound chemistry face states with a scaling limiter (NaN at finite-cylinder rims) - #1975

Contribution Policy

All contributions are expected to demonstrate:

• A clear understanding of the codebase
• Alignment with product direction
• Thoughtful reasoning behind changes
• Evidence of real-world usage or hands-on experience with the problem

If these expectations are not met, we would prefer to implement the changes ourselves rather than spend time reviewing low-effort submissions.

Fixes #1964

Root cause

This is not a levelset or image-point bug. I checked every ghost point of the repro cylinder offline (same marker, levelset, image-point and stencil logic as master). Every image point lands in fluid, and every interpolation stencil has at least 4 fluid points (cm3c_cases/debug_E/cyl_rim/tools/geom.py).

The failure happens like this:

1. At the cap/side rim, the ghost column under the first side-fluid cell (axial depth 0.5 dx) mirrors different faces. Layer 1 (radial depth 0.16 dx) is nearest the side and mirrors the expansion: p = 1.27e5. Layers 2–4 are nearest the cap and mirror the stagnation region in front of it: p = 2.3–2.8e5.
2. With a hot isothermal wall (thermal_bc = 1, Twall = 1800 K), the ghost densities in that column are 0.14, 0.24, 0.27 and 0.28 kg/m^3, next to fluid at 1.05. The column is non-monotone.
3. WENO picks the smooth all-ghost sub-stencil and extrapolates the inward-rising ghost density out to the fluid/ghost face. The face density undershoots to 0.059, so face T = p·W/(ρR) ≈ 7600 K.
4. Above the NASA-7 fit (3500 K for O2/CO2 in this mechanism), cp_O2/R evaluates to −10. HLLC then gets a NaN sound speed and returns a NaN flux on that face. Only the fluid cell goes NaN, because IB-cell RHS is zeroed.

Instrumented trace from a restart at step 21 of the repro: at stage 3 of step 22, face (12, 17.5, 24) has ρ_R = 0.059 and p_R = 1.31e5, and the y-flux there is NaN. Non-reacting finite rods (no chemistry, adiabatic wall) run cleanly. With chemistry, an adiabatic IB (thermal_bc = 0) or an ambient-temperature wall (Twall = 350 K) also runs cleanly. The hot isothermal wall is the ingredient that drives this to NaN.

A convex corner has no single ghost value that suits both faces. A real fix at the IB level needs a new corner treatment (e.g., per-direction ghost states). This PR bounds the reconstructed face states where the NaN is produced.

#1969 alone does not fix this. #1969 extrapolates the hot-wall ghost from where the image point is actually sampled. On the repro it cools the shallowest rim ghost layer (step 21: 3176 K → 2063 K, ρ 0.137 → 0.193) but leaves layers 2–4 at ~3120 K, so the column is still non-monotone. The NaN is delayed, not removed: the 3D repro NaNs at step 26 instead of 22 (cell (12, 18, 21), the same rim position), and the 2D rectangle case at step 70 instead of 24.

Change

• m_rhs.fpp: new s_bound_face_states(id), run after reconstruction when chemistry is on. It is a Zhang–Shu scaling limiter. Each cell's two face states (all primitive variables) are scaled toward the cell state by one θ ∈ [0, 1]. θ is the largest value that keeps face ρ, p and R_mix above 1e-12 of the cell value and face T = p/(ρ R_mix) ≤ T_surface_max. T_surface_max is 5000 K, the thermo window m_ibm already uses. ρ, p and R_mix are all linear in θ. So the positivity bound is closed form, and the T bound is the first root of a quadratic. That makes θ continuous in the data, with no iteration. θ = 1 (no change) wherever the faces are already in bounds. There is no lower T bound: an earlier 200 K bound moved 1D -> Chemistry -> MultiComponent Diffusion. The range and the GPU loop style follow the species-bounds limiter (Keep species mass fractions in [0, 1] by construction instead of clipping #1963), so it runs over the range the reconstruction filled, with θ shared by both faces of a cell.
• This replaces the previous version of this PR, which reverted an offending face to first order. That was a hard switch on one face; this is bounded, continuous and keeps the reconstruction wherever it is admissible.
• New test 2D -> Chemistry -> IBM -> Hot Wall Rectangle Corners (4646053A). It uses h2o2.yaml, 72x48 cells, 30 steps, ~3 s on CPU. It is the 2D analog: Mach 1.5 air onto an 1800 K isothermal IB rectangle. On master it aborts at step 30 (ICFL is NaN).

Evidence (CPU, CCE, serial --no-mpi)

• 3D repro (x-aligned finite cylinder, D = 2e-4, L = 3D, 12 cells/D, 72x48x48, Mach 1.5 axial air, viscous, no-slip, chemistry with reactions off, surface reactions off, Twall = 1800 K, fixed dt = 5.31e-9; run in two restart segments of 150 steps):
  • RED on master: NaN at step 22 at cell (12, 17, 24), the first fluid cell outside the side, 0.5 dx behind the upstream cap.
  • GREEN with this PR: all 300 steps, exit 0, no NaN. The limiter acts on 154 distinct cells (144 IB ghost cells, 10 fluid cells), all at the rims. Minimum θ is 0.11; unlimited face T reached 7e6 K in ghost cells and 1.5e4 K in fluid cells. Near the rim, the hottest fluid cell at step 300 is 2553 K. The previous first-order revert gave 2852 K there.
• 2D rectangle, 300 steps (same setup as the test): master NaNs at step 24; this PR runs all 300 steps. The limiter acts on 20 cells, all at the corners (19 IB ghost cells, 1 fluid cell); minimum θ is 0.27. In the 30-step test it acts on 14 ghost cells only.
• ./mfc.sh test --no-mpi --no-gpu:
  • -o IBM: 62 passed. A69D2D28 failed once with a transient No such file or directory: '0' and passed on rerun.
  • -o Cylinder: 8 passed.
  • -o Chemistry: 23 passed.
• Goldens regenerated:
  • 4646053A: the new test. Its golden was made with the revert and differs from this limiter's result.
  • B317404C (2D -> IBM -> Vieille Burn Rate) and E5B66084 (2D -> Example -> ibm_burning_grain): the same case. On master, a blowing ghost cell inside the grain reconstructs a negative face density (ρ_face = −0.059 next to a cell ρ = 0.020, so face T < 0). HLLC is fed that state on master. The revert only touched that face, and the change stayed under the 1e-3 test tolerance. The limiter scales both faces of that cell (θ = 0.20). Fluid momentum moves by up to 3.7% relative where |value| > 1e-6, and the IB force (small in magnitude) changes a lot. With a per-face θ instead, I checked that B317404C stays within its tolerance. I kept the shared per-cell θ for consistency with Keep species mass fractions in [0, 1] by construction instead of clipping #1963. Reviewers may prefer per-face θ.
  • No other golden moved.
• ./mfc.sh precheck passes.

Not verified

• GPU builds (OpenACC/OpenMP) and MPI runs of the new kernel.
• Compilers other than CCE. I set the new test's tolerance to 1e-5 instead of 1e-10: its corner faces amplify roundoff, and CI lanes for Keep species mass fractions in [0, 1] by construction instead of clipping #1963's similar test differed from the CCE golden by up to 6e-7 relative.
• The original reacting-surfaces rod case (surface reactions on). That is on a fork branch, not master.
• Flows that are legitimately above 5000 K get their face states scaled toward the cell state. Cell states are not limited.
• Rim fluid still heats to ~2500 K, above Twall. That is the corner ghost treatment, which this PR does not change.